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Abstract. Numerical algebraic geometry uses numerical data to de- 
scribe algebraic varieties. It is based on numerical polynomial homotopy 
continuation, which is a technique alternative to the classical symbolic 
approaches of computational algebraic geometry. We present a pack- 
age, whose primary purpose is to interlink the existing symbolic meth- 
ods of Macaulay2 and the powerful engine of numerical approximate 
computations. The core procedures of the package exhibit performance 
competitive with the other homotopy continuation software. 

Numerical algebraic geometry |15|. [16] is a relatively young subarea of 
computational algebraic geometry that originated as a blend of the well- 
understood apparatus of classical algebraic geometry over the field of com- 
plex numbers and numerical polynomial homotopy continuation methods. 
Recently steps have been made to extend the reach of numerical algorithms 
making it possible not only for complex algebraic varieties, but also for 
schemes, to be represented numerically. What we present here is a descrip- 
tion of "stage one" of a comprehensive project that will make the machinery 
of numerical algebraic geometry available to the wide community of users of 
Macaulay2 [9], a dominantly symbolic computer algebra system. Our open- 
source package dubbed NAG4M2 [11 J and NumericalAlgebraicGeometry [9j 
was first released in Macaulay2 distribution version 1.3.1. 

"Stage one" has been limited to implementation of algorithms that solve 
the most basic problem, upon solution of which the majority of other prob- 
lems depend: 

Given polynomials fi, ■ ■ ■ , f n £ C[x±, . . . , x n ] such that they 
generate a O-dimensional ideal I = (/i, . . . , f n ) find numer- 
ical approximations of all points of the underlying variety 
V(I) = {x | /(x) = 0}. 

This task is accomplished by applying the idea of homotopy continuation. 
To solve a target polynomial system / = • • • , f n ) = construct a start 
polynomial system g = (gx, . . . , g n ) with a "similar structure" (the meaning 
of this will be explained later), but readily available solutions. Define a 
homotopy, 

(1) h=(l-t)g + jtfeC[x,t], 7€C*, 

which specialized to the values of t in the real line interval [0, 1] provides a 
collection of continuation paths leading from the (known) solutions of the 
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start system, g = h\ t=0 , to the (unknown) solutions of the target system, 
/ H ,• 

One option for a start system with the "similar structure" and readily 
available (regular) solutions, g = (x^^ 1 — 1, . . . , x^ g ^ n — 1), leads to the 
so-called total degree homotopy, for which the following statement enables 
numerical computation. 

Theorem. For all but finitely many values of 7 in the homotopy |7]) the 
homotopy continuation paths have no singularities with a possible exception 
of the endpoints corresponding to t = 1 . 

Every solution of the target system, provided there are finitely many, is 
an endpoint (at t = 1) of some continuation path. 

Proof. Let P be the linear space of systems of n polynomials in n variables 
with complex coefficients. Note that P = C m for some m E N and systems 
with at least one singular solution form a Zariski closed set A C P. With 
the homotopy ([I]) we satisfy the conditions of Lemma 7.1.3 of [16] with an 
exception of the direction of the homotopy (in |16] the parameter t varies 
from 1 to 0). The Lemma concludes that for all but finitely many choices 
of 7 the homotopy ([I]) for t G [0, 1) misses the set A. □ 

The homotopy continuation argument above belongs to the core of clas- 
sical algebraic geometry; it has been known in the beginning of the 20th 
century. 

Differentiating the homotopy equation h = gives the following system 
of ODEs 

<2) 

where h x is the Jacobian of the homotopy (with respect to x) and ht is 
the derivative with respect to the parameter t. The solutions of ([2]) for 
t 6 [0,1] with initial conditions given by the solutions of the start system 
are the continuation paths we need. Finding continuation paths approxi- 
mately, therefore, reduces to numerical solving of systems of ODEs. There 
is an additional advantage: at any point t = to we can refine an approxi- 
mate solution to the polynomial system h\ t=tQ = with Newton's method. 
Provided the numerical tracking procedure has not deviated from the given 
path, this brings the approximation as close to the path as desired. 

The basic tracker operates by alternating predictor (a numerical integra- 
tion step) and corrector (several applications of Newton's operator) as shown 
in Figure [T] The ultimate goal of the tracker is to approximate the end of a 
continuation path given an approximation of its beginning. 

An introductory example of tracking a total-degree homotopy is below: 



il 

i2 
i3 
i5 



needsPackage "NumericalAlgebraicGeometry" 
R = CC[x,y] ; 

S = {x~2-l,y~2-l}; T = {x~2+(y-5)~2-16, x*y}; 
solsS = {(1,-1) ,(1,1), (-1,1), (-1,-1)}; 
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16 : track(S ,T, solsS) 

06 = { [M,t=. 04762] , {2.955e-17, 1}, {9.912e-14, 9}, [M,t= . 04762] } 

17 : track(S,T,solsS,gamma=>0.6+0.8*ii) 

07 = {{1.374e-13-4.782e-14*ii, 9}, {3*ii, -5.727e-17+9.943e-17*ii}, 

{-2.026e-17-2.601e-17*ii, 1}, {-3*ii, 4 . 792e-18-3 . 185e-19*ii}} 

Note that the functions in the example use 7 = 1, the default value of 
the parameter 7 in ([!]) , which often results in a homotopy that goes through 
a singular point when both start and target systems have real coefficients. 
In the example above the tracking of two paths fails as ^ 1^0 04762 = 
has singular solutions. Picking a random complex value for 7 (in practice, a 
number is chosen on the unit circle with the uniform probability distribution) 
results in a regular homotopy with probability 1 according to the Theorem. 
In fact, the black-box solver solveSystem does exactly that when it calls 
track. 

Software implementation strategy. Despite the seeming simplicity of 
the approach described above, there are many technical details that an im- 
plementer has to think through: choosing a predictor mechanism, dynamic 
step adjustment, etc. The development of NAG4M2 has started with de- 
veloping the basic function called track that given a target system, a start 
system, and a list of start solutions numerically tracks the corresponding 
homotopy paths as the continuation parameter t is varied from to 1. This 
has been carried out in the top level language of Macaulay2; more than a 
dozen optional parameters supplied to track allow users to experiment with 
various settings of the tracker and set numerical thresholds that control the 
precision of the computation. 

The language of Macaulay2 is, however, an interpreted language; the per- 
formance of the code was far from optimal. That is why it was crucial to 
implement the computationally intensive parts of the code in C++ in the 
Macaulay2 engine. At the present moment pieces of the source code are 
written in three languages: 
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• Macaulay2 language: processing the input and setting up evalua- 
tion, predictor, and corrector routines; launching the routines in the 
engine; managing the output; 

• C++ language: fast execution of the predictor-corrector steps using 
standard double floating point arithmetic; 

• D language: providing interface, converting objects of Macaulay2 
top level classes into objects of types and classes used by the C++ 
engine. 

The C++ part relies on lapack pQ for numerical linear algebra. This is the 
only external library that has been used so far. 

Performance. In addition to native implementation the user is given an 
option to outsource the computation to one of the three polynomial homo- 
topy continuation softwares: Bertini [2], Hom4PS2 [10], or PHCpack [17j . 
This depends on the availability of the software for the user's platform. 

One of the goals of "stage one" of the project was to achieve performance 
competitive with the existing software. The following timings (obtained on 
a single core of a 64-bit Linux system) demonstrate the results for several 
test problems with a moderate number of solutions (less than 10,000): 

• Random^: a system of dense polynomials of degree d in n variables 
with random coefficients; 

• Katsura n : a classical benchmark with one linear and n — 1 quadratic 
equations in n variables; 

• GEVP n : the system corresponding to a generalized eigenvalue prob- 
lem, Av = XBv for n x n (randomly generated) matrices A and 
B. 

All these problems have regular solutions and do not encounter near- 
singularities, i.e., the continuation paths are sufficiently far from each other, 
so that double-precision arithmetic is enough to carry out the computations. 
All runs for all software are made in standard double precision with default 
settings rj 
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For all systems except the last one, the number of solutions equals the 
total degree, while for GEVP n we supply an optimal homotopy with a start 



Disclaimer: The timing ranks may depend on our selection of benchmarks. In fact, there 
are many problems which can not be solved with one program, but are solved by another 
quickly. In addition, timings depend on a variety of factors including the hardware, the 
operating system, and the parameters of continuation. 
2 HomJ h PS2 has no option to handle a user-defined homotopy. 
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system having exactly n solutions. One may rerun the examples used to 
obtain the timings using the files showcase. m2 and benchmarks . m2 in the 
Macaulay2 repositorjj^j The (current) default options of the function track 
were used to complete the tests. In particular, Predictor => RungeKutta4 
specified a predictor of the fourth order, which performed slightly better than 
lower order numerical integration algorithms on all mentioned examples. 

Technical details. There are many factors that may affect performance of 
the homotopy tracking procedure: the choice of predictor, parameters that 
control dynamical step adjustment, efficiency of linear algebra subroutines, 
etc. 

One common bottleneck in the computation could be the evaluation of 
polynomials. All test examples mentioned above, except Katsura n , are 
evaluation-intense. The evaluation of polynomials presented in the dense 
form, which is the form Macaulay2 uses, is quite expensive. To speed up 
the computation we employ an extended Horner scheme (see e.g. [ZJ), a 
generalization of the Horner scheme for evaluating univariate polynomials. 
While not necessarily an optimal scheme (in the sense the univariate Horner 
scheme is) it makes an attractive design choice for two reasons. First, we can 
encode the evaluation procedure in a straight-line program (SLP), a program 
evaluating which can be achieved without branching or looping. An SLP 
corresponding to an extended Horner form is often much shorter than that 
for the dense form. Second, the automatic differentiation of polynomials 
represented in Horner form is convenient as the size of an SLP representing 
the evaluation of a polynomial function and its derivatives simultaneously 
is not much larger than the size of an SLP for the polynomial function itself 
(see the corresponding theoretical complexity result in [3]). 

We experimented with speeding up the evaluation of SLPs even further 
by compiling SLPs at runtime. This results, on some evaluation-intense 
examples, in up to 10 times faster execution]^] At the moment, the option of 
compiled SLPs is not developed for all platforms and needs to be examined 
further. For example, the current compilation time may exceed the gain in 
computation time in some cases. 

Certification. One important feature that distinguishes NAG4-M2 from 
other software is that of the certified homotopy tracking procedure devel- 
oped in joint work with Beltran [5]. In general, the heuristic algorithms dis- 
cussed above need careful tuning before path-jumping (landing on a wrong 
path) is eliminated. The most commonly used step control procedures are 
described in [3] along with Bertinis adaptive multiprecision approach for 
overcoming poor conditioning in numerical linear algebra while tracking the 
near-singular continuation paths. 



svn : //macaulay2 .math .uiuc . edu/Macaulay2/trunk/M2/Macaulay2/packages/NumericalAlgebraicGeometry 

4 This trick seems to be used in Hom4PS2 as well. However, the timings reported for the 
test examples above are obtained by running all programs without runtime compilation. 
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In contrast, a certified homotopy tracking algorithm makes a safe choice 
for a size of step in the continuation, guaranteeing that an approximate 
zero created in the next step is associated to the same homotopy path. If all 
computations are carried out in exact arithmetic the certified algorithm gives 
a rigorous proof of the results obtained. Our implementation of the certified 
tracking [3] is invoked by passing the option Predictor=>Certif ied to 
the function track. At the moment the linear algebra computations are 
carried out in the standard double precision, thus stopping one step short 
of rigorous certification. The work accomplishing this last step using exact 
linear algebra and the robust a-theory is underway [6]. 

Future. The package in its current state provides the base for the higher- 
level numerical algebraic geometry routines. Amongst them are robust com- 
puting of approximations to singular solutions using endgames [TBI §10.3] 
and deflation [131 E] ■> irreducible decomposition of positive-dimensional va- 
rieties |16} §15], numerical primary decomposition |12| . etc. 

Homotopy tracking with arbitrary precision, which both Bertini and PHCpack 
implement to various extents, should be available in Macaulay2 in the fu- 
ture. Arbitrary precision floating point arithmetic is already in place in 
Macaulay2 (via the MPFR library [8]); however, a fast linear algebra im- 
plementation on the level of the engine is necessary for solution refinement 
and tracking near-singular paths. 

Due to independence of tracking procedures for any collection of homo- 
topy paths most of the algorithms in numerical algebraic geometry scale well 
if parallelized. Since the amount of data that needs to be communicated be- 
tween CPUs is small in relation to the computational costs, the tasks can be 
distributed over heterogeneous clusters with slow interconnect. For the large 
problems these properties give homotopy continuation an edge over Grobner 
bases techniques, which are currently the main engine of Macaulay2. 

We will explore ways to convert numerical results into symbolic results and 
vice versa. For instance, the algorithms in the package will be able to sample 
points on a component of a variety with an arbitrarily high precision, which 
could be used to construct the defining ideal of the component bypassing 
the need for (symbolic) primary decomposition. 
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